function [res C] = ptrSort (varargin)

   ## usage:  [res C] = ptrSort (varargin)
   ##
   ## sort x according to correlation to Y.z
   ## opt indicates single or grouped (acc. to variables V) operation

   x = varargin{1} ;
   J = varargin{2} ;
   y = varargin{3} ;
   if (length(varargin) < 4)
      opt = "single" ;
   else
      opt = varargin{4} ;
   endif 

   if isfield(y, "isP")
      J = y.isP ;
      res1 = res2 = C1 = C2 = [] ;
      if any(J)
	 varargin{3} = mean(y.z(:,J), 2) ;
	 [res1 C1] = ptrSort(varargin{:}) ;
      endif 
      if any(!J)
	 varargin{3} = mean(y.z(:,!J), 2) ;
	 [res2 C2]= ptrSort(varargin{:}) ;
      endif 
      res = merge_perm(res1, res2) ; C = [C1(res1); C2(res2)](:) ;
      return ;
   endif 

   switch opt

      case {"single"}
 
	 C = abs(corrcoef(y, x)) ;
	 # C = re(zscore(x), zscore(y)) ;
	 [s res] = sort(C, "descend") ;

      case {"multi"}

	 Jv = values(J)' ;

	 jj = 0 ;
	 for j=Jv
	    [b, bint, r, rint, stats] = regress (y, x(:,J==j)) ;
	    S(++jj,:) = stats ;
	 endfor 
	 [s J1] = sort(S(:,1), 1, "ascend") ;
	 [s J2] = sort(S(:,2), 1, "ascend") ;
	 [s J3] = sort(S(:,3), 1, "descend") ;
	 [s J4] = sort(S(:,4), 1, "descend") ;

	 J1(J1) = J2(J2) = J3(J3) = J4(J4) = (1:length(Jv)) ;
	 [s p] = sort(mean([J1 J2 J4], 2), "descend") ;

	 res = C = [] ;
	 jj = 0 ;
	 for j=Jv
	    res = [res find(J == Jv(p(++jj)))] ;
	 endfor 

      otherwise

	 fid = fopen("~/tmp/log", "w") ;
	 p = C = [] ; r = y ;
	 for j=1:150
	    [r, p] = rrgr(r, x, p) ;
	    fmt = [repmat("%7d", 1, length(p)) "\n"] ;
	    fprintf(fid, fmt, p) ;
	    fflush(fid) ;
	 endfor 
	 fclose(fid) ;
	 res = p ;

   endswitch 

   if ((nu = sum(isnan(C))) > 0)
      res = [res(nu+1:end) res(1:nu)] ;
   endif

endfunction


function I = merge_perm (I1, I2)

   ## usage:  I = merge_perm (I1, I2)
   ##
   ## merge two permutations

   if (isempty(I1) || isempty(I2))
      I = [I1 I2] ;
      return ;
   endif 

   n = length(I1) ;
   wI = 1:n ;

   w(2*wI-1) = I1 ; w(2*wI) = I2 ;

   I = [] ;
   while !isempty(w)
      I = [I w(1)] ;
      w = w(w != w(1)) ;
   endwhile 

endfunction
